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I.   INTRODUCTION 

Two  basic  methods  exist  for  solving  the  response  of  a 
system  which  generally  is  presented  as  an  integro-differential 
equation:   analytic  methods  and  numerical  methods. 

For  simple  problems,  analytic  methods  give  precise 
solutions  and  have  the  advantage  that  these  solutions  can 
be  given  as  a  function  of  the  parameters  of  the  system. 

As  the  order  of  the  differential  equation  increases, 
numerical  methods  give  less  labor  in  providing  a  solution 
and  can  be,  at  some  extent  the  only  method  available. 

Laplace  transform  is  limited  to  the  case  of  linear 
differential  equations.   Nonlinear  differential  equations 
are  out  of  the  scope  of  this  technique,  preventing  its 
application. 

EVen  sometimes  with  linear  systems,  if  the  inverse 
integral  is  complex,  it  is  difficult  and  sometimes  impossible 
to  find  analytically  the  time  function. 

However,  it  is  possible  to  approximate  the  Laplace 
transform  of  a  linear  system  by  a  z-transform,  obtaining 
thereafter  a  time  series  representing  the  values  of  the 
response  at  equally  spaced  intervals  of  time,  by  synthetic 
division. 

The  method  [1]  described  and  used  In  this  paper  utilizes 
the  z-transform  directly  in  order  to  obtain  the  solution  of 
the  linear  and  nonlinear  systems. 


In  this  paper,  an  explanation  is  given  of  this  method 
applied  to  linear  systems  and  its  extension  to  time-varying 
and  nonlinear  systems. 

A  computer  was  used  for  faster  and  more  reliable  results, 
but  the  method  described  is  compatible  with  handwork. 

C.P.U.  time  in  the  computer,  for  this  technique,  is 
much  smaller  than  other  numerical  methods  due  to  the  less 
number  of  iteration  steps  needed.   This  technique  could 
reduce  solution  times  by  several  orders  of  magnitude  over 
conventional  numerical  methods. 


II.   PRESENTATION  OF  THE  METHOD 

A.   PROCEDURE 

The  basic  procedures  of  the  method  are  given,  before 
the  description  or  the  theory,  to  show  its  simplicity. 

Suppose  that  one  has  the  Laplace  transform  of  an  output 
response.   It  will  be  usually  in  the  following  form: 


a  +  a,  s  +  a0s2  +  •••  +  a,  sk  +  ...  +  a  sn 

PCs)  =  -° 1 2—2 *-£ 2_      (1) 

b  +b,s+b~s^  +  •  •  •  +  b,  sK  +  ...  +  b  s 
o   1    2  k  m 


The  steps  to  obtain  the  time  response  are: 

1.  Express  the  function  F(s)  as  a  rational  fraction  in 
powers  of  s~  by  dividing  the  numerator  and  denominator 
by  sm. 

2.  Substitute  for  s~   a  rational  fraction  in  powers  of  z~ 
obtained  from  Table  I  and  rearrange  F(s)  as  a  rational 
fraction  in  powers  of  z   . 

3.  Divide  the  resulting  expression  by  T,  where  T  is  the 
time  interval  between  points  at  which  the  result  is  desired. 
^ .   Expand  the  fraction  by  synthetic  division  into  a  series 
of  the  form: 

DQ + D1z~1 + D2z~2  +  ...  +  Dnz~n  +  ...         (2) 


where  D^ ,  the  coefficient  of  z   ,  is  the  approximate  value 


of  the  time  response  at  t  =  nT. 
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B.   LINEAR  SYSTEMS 

Take  an  arbitrary  function  f(t)  as  in  Fig.  1(a).   Its 
Laplace  transform  and  inverse  Laplace  transform  are  defined 
as : 


F(s)  =  /  f(t)  e"st  dt 


(3) 


f(t)  = 


C+j' 


St 


2ttj 


/    f(s)  eou  ds 


c-j< 


(4) 


y  ■ 


|&)- 


s  - 


i 
i 

0 


/w*a.=>|M 


T  2.T  3? 


"i 1 1 r 

*l7 


?  ZP  3' 


Fig.  1(a):   Arbitrary  f(t)      Fig.  1(b):   Sampled  f(t) 


Figure  1(b)  shows  f(t)  sampled  by  shifted  unit  impulses, 
The  result  is  a  series  of  impulses  at  t  ~  nT. 


Applying  Laplace  to  the  sampled  function  one  gets 


F*(est)  =  Z        f(nt)  e~snT  (5) 

n=0 

sT 
substituting  in  (5)  e   by  z  so: 


z  =  e  (6) 


the  result  obtained  is  known  as  the  z-transform: 


F*(z)  =   Z    f(nT)  z  n  (7) 

n=0 


and  the  inverse  is  given  by: 

f'(nT)  =  £   P*(z)  zn-1  dz  (8) 

where  the  contour  is  a  circle  centered  at  the  origin  and 
including  all  the  singularities  of  F*(z). 

The  equation  (6)  maps  the  j to  axis  of  the  s-plane  into  a 
unit  circle  around  the  origin  in  the  z-plane. 

There  are  other  methods  for  obtaining  direct  and  inverse 
z-transforms. 

In  closed  form  the  function  will  appear  as: 


a  +  a,  z~  +  a0z~  +  •••  +  a  z  n 

"<«)  -   "    1    -!    2    -2 —  (5> 

bn  +  bn  z  L  +  b0z   +  • ••  +  b  z 

0    12  m 
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Division  of  the  numerator  by  the  denominator,  will 
result  in  the  coefficients  of  f(nT)  as  defined  in  (7), 
without  the  use  of  the  inversion  integral  formula  (8). 

In  order  to  apply  this  useful  characteristic  of 
z-transform  theory  to  non-sampled  systems  one  uses  the 
inverse  Laplace  integral  (4), 


c+ j00 
f(t)  =  2YT    /  '   PCs)  est  ds 


which  can  be  expressed  piecewise  as: 


,         Jtt/T  J» 

f(t)    -  JU-        /  f(s)eStds   +   i-        /      [P(s)eSt  +f(-s)e"St]ds 

^J    -Jtt/T  ^J    Jtt/T 


(10) 


The  way  the  contour  is  divided  is  shown  in  Fig.  2.   The 
first  integral  gives  the  part  of  the  contour  shown  between 
-(tt/T)  <_  03  <   tt/T  .   The  remaining  contour  w  >  (tt/T)  , go  <  -(tt/T) 
is  given  by  the  second  integral. 

If  1/T  is  sufficiently  large,  a  good  approximation  will 
be  held  by  the  first  integral  of  (10)  and  the  second 
integral  may  be  discarded,  being  the  resulting  error  given 
by: 


J°° 
/ 

-Jtt/T 


e(t)  =  2t7J     ;    f(s)eSt  +  F(-s)e-st  ds       (11) 
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f 


-jyt 


c< 


Figure    2:      Contour   in   s-plane   for   approximate    solution 


The   approximation,    f.(t),    of   f(t)    is   given  by 


fA<*>    =    2¥J 


Jtt/T 


st 


/  F(s)    eou   ds 


-Jtt/T 


(12) 


Making  t   =   nT 


fA(nt)    =    2t7J 


jTr/T  snT 

/  P(s)    eSnl    ds 

-Jtt/T 


(13) 


Solving  eq .  (6)  for  s,  one  obtains 


s  =  1/T  In  z 


(14) 
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Making  the  substitution  of  (14)  into  (13)  the  following 
equation  is  obtained 


VnT)  ■  2?J   J7/T  PC1/T  in  z)  enT(l/T)ln  z  d(l/T  In  Z) 
-Jtt/T 

(15) 


fA(nT)  =  -jfj  <*  |  P(l/T  In  z)  zn_1  dz  (16) 


where  the  contour,  again,  is  the  unit  circle  in  the  z-plane. 

If  eq .  (16)  is  compared  with  eq.  (8),  it  is  seen  that  it 
has  the  form  of  an  inverse  z-transform  except  for  the  constant 
multiplier  1/T. 

On  the  other  hand,  as  P(l/T  In  z)  is  a  transcendental 
function,  it  is  not  expandable  directly  into  a  series  of 
powers  of  z~  by  synthetic  division.   The  synthetic  division 
technique  can,  however,  be  used  if  a  proper  approximation 
could  be  found. 

If  F(s)  is  expanded  in  a  series  of  descending  powers  of 
s  like 

F(s)  =  a1s~1  +  a2s~2  +  •••  (17) 

a  solution  of  the  form 

t2 
f(t)  =  a1  +  a2t  +  a  |y  +  •••  (18) 

can  be  obtained,  but  this  solution  is  not  always  valid. 
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The  way  to  follow  Is  to  rearrange  F  (l/T)/ln  z  as  a 
rational  function  of  two  polynomials  In  powers  of  s   =  T/ln  z. 
For  this  purpose,  the  power  series  for  In  z  used  is: 

In  z  =  2[u  +  1/3-u3  +  1/5-u5  +  •••]  (19) 

where : 

1  -  z-1 
u  =  


1+  z"1 


So  the  following  equation  is  obtained; 


1  =   T   =  T/2 

s   In  z     ,  i  /0   3  ,  n/r   5 
u  +  1/3  uJ  +  1/5  u 


(20) 


By   synthetic   division   the    following  Laurent   series   is 
obtained: 


1  _       T  T   r   1        1  ^       3         ^^       5    ,  n  /o-n 

b      WT  =  2  [  u  -  3  u  "  TH5"  u        P5  u    +]  (21) 


To  obtain  series  expansion  of  s   ,  bothsides  of  (21) 
can  be  raised  to  the  k   power. 

Table  I  (on  page  15)  was  obtained  by  using  only  the 
principal  part  and  the  constant  term  of  the  Laurent  series 
The  result  is  of  the  form 


1H 


TABLE    I 


1  —  ■"■    - 

s-k 

z-Form  =  Fk(z  1) 

s"1 

T  1  +  z"1 
2  1-z-1 

s-2 

2        -1    -2 
T^  1  +  lOz  ±  +  z 

12    (1-z"1)2 

s-3 

T3  z"1  +  z"2 
2  (l-z-1)^ 

a"1* 

T"  z"1  ♦  Hz"2  ♦  z"3    t" 

s 

6    (1.B-l,t      YzO 

s-5 

T5  z"1  +  llz-2  +  llz~3  +  z 

2*                  (l-z-1)^ 
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(1  -  z   ) 


(22) 


where  N.  is  a  polynomial  in  powers  of  z   .   The  right  hand 
side  of  (24)  is  called  the  "z-form"  for  s~k. 


Example  1  -  Third  Order  System  (Fig.  3) 


ft&J 


± 


o- 


i/scs+i; 


a,-- 


0    s3+s-+s 


Fig.  3:   Typical  Feedback  System  for  Example  1 

With  a  step  input,  as  shown,  the  Laplace  transform  of 
the  output  is : 


VS)  =  -J^~2 

u      s^  +  s  +  s 


(23) 


Following  the  step  by  step  procedure  as  given  in  the 
Introduction, 


Vs)  =  — h — ^ 

u      1  +  s  x  +  s  J 


(24) 
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— k 
Now,  substituting  the  correspondent  forms  of  s    from 

Table  I  and  dividing  the  result  by  T: 

6  *(z)  =  6T2(z-l  +   z-2) ^ 

0A      (12+6T+T2)  -  (36+6T-9T2)z_1  +  (36-6T-9T2)z~2  -  (12-6T+T2)z~3 

(25) 

The  solution  is  obtained  by  choosing  T,  the  time  interval 
at  which  rate  our  solution  is  obtained,  and  proceeding  in 
dividing  the  denominator  into  the  numerator. 

Let  us  assume  T  =  1.0,  resulting  in  the  following  output 


60A*(z)  -  6Z"'  +  6Z"!2 rs  (26) 

u        19  -  33z  x  +  21z  d  -  7z  5 


By  long  division,  the  result  is: 

.3l6z_1 +  ,864z"2  +  1.15z~3  +  l.l6z~  +  1.06z"5  +  .987z~D  +  •  •• 

The  accuracy  improves  if  smaller  value  for  T  is  chosen. 
Letting  T  =  0.5  in  this  same  example,  substituting  on  eq.  (25) 

6  *(z)  =  1.5^  +  1.5z"2  (27) 

UA       15.25  -  36.75z-1  +  30.75z  d  -   9.25z  3 

carrying  out  the  long  division,  the  result  is: 

0.984z_1  +  .335z~2  +  .  6lOz"3  +  .  853z_i|  + 

The  points  obtained  are  plotted  together  with  the  precise 
solution  in  Fig.  4. 
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Figure  4 


X-scale 
Y-scale 


2.0  units/inch 
.5  units/inch 


—  Precise  solution 
x   Approximate  solution 

for  At  =  1.0 
+   Approximate  solution 

for  At  =  0.5 


Example  2  —  Solution  of  Linear  Differential  Equation 

To  use  the  present  technique  in  obtaining  the  solution 
of  integro-differential  equations,  one  takes  the  Laplace 
transform  of  the  equation,  initial  conditions  included  as 
usual,  and  then  proceeds  as  before.   However,  after  trans- 
forming the  equation,  the  following  rules  can  be  applied: 
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1.  If  the  degree  of  the  denominator  of  the  Laplace 
transform  is  higher  than  the  degree  of  the  numerator  by  two 
or  more,  proceed  directly. 

2.  If  the  degree  of  the  numerator  is  one  less,  equal  to, 
or  greater  than  the  degree  of  the  denominator,  divide  until 
the  remainder  meets  the  criterion  of  Rule  1. 

Using  the  procedure  stated  in  the  rules  above  one 
eliminates  the  high  frequency  from  P(s),  that  is  implied  if 
the  numerator  is  one  degree  less  than  the  denominator. 

To  illustrate  this,  a  simple  example  is  shown: 


f|  +  Y  =  1   .  (28) 


Assuming  for  initial  condition  Y(0)  =  u,  the  Laplace 
transform  Y(s)  is: 


*<s>  ■  Ix^iy  -  7777  <2*) 


Using  Rule   2    above,    one   gets 


1  -u 


Y(s)    -  g+      X"U  (30) 

s    +  s 


In  this  procedure,  the  step  function  of  amplitude  u  is 
added  to  the  result  obtained  after  operating  upon 
(1  -u)/(s2  +  s). 
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Now,  letting  u  =  0.1,  and  expressing  the  second  term 
of  eq.  (30)  in  powers  of  s   , 


-2 


yl(8)  =  0^-_ 


(3D 


1  +  s 


Substituting  the  z-forms  of  Table  I  (on  page  15)  and 
dividing  by  T, 


v  */^  _     0.9T(1  +  lOz"1  +  z  2) 

XK       (12  +  6T)  -  24z  x  +  (l2-6T)z 


(32) 


Letting  T  =  0.1  and  carrying  out  the  long  division  one 


gets 


,00714  +  .085z  -1-  +  .l63z  2  +  .233z~3  +  .296z-i4  + 


(33) 


This   result    is   plotted   in   Fig.    5   after   adding  the    step 
function   of   amplitude    0.1. 

millllllllllllllllllllllllllllllllllllllllllllllllllll 
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—  Precise  solution 
+  Approximate  solution 
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C.   TIME-VARYING  SYSTEMS 

A  linear  differential  equation  with  time-varying 
coefficients  may  be  expressed  as 


3n 


dt 


n 


fn(t)Y 


,n-l 


dt 


n-l 


fn-l(t)Y 


+  f0(t)Y  =  f(t)    (3*0 


where  f  (t),  f  , (t) ,. . . ,fQ(t)  are  functions  of  the 
independent  variable  t. 

To  solve  this  type  of  equation  the  following  steps  may 
be  executed: 

1.  Obtain  the  Laplace  transform  of  (3*0  in  the  usual 
manner  including  initial  conditions. 

2.  Proceed  as  in  the  constant  coefficient  case  discussed 
above.   In  performing  the  long  division  process,  however 3 
change  the  values  of  C  ,  ...,Cn  at  each  step  in  the  division 
process  to  the  values  at  the  corresponding  times  of  the 
functions 


Cn  =  fn(t) 


CQ  =  fQ(t) 


(35) 


In  general 


c1  =  f±(t) 


(36) 
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D.   NONLINEAR  SYSTEMS 

A  similar  manner  like  used  in  II. C.  may  be  applied  to 

nonlinear  differential  equations.   The  regression  equation 

is  obtained  by  long  division  process  as  in  the  time  varying 

example. 

However,  the  factors  C  ,  C   ,,...,  Cn  must  now 
3  n'  n-1     '   U 

represent  functions  of  the  dependent  variable. 

Generally,  a  nonlinear  differential  equation  can  be 
represented  as: 


,n 


dt 


n 


f  (Y  —     )Y 
Jnv  »dt»" '• ; 


an-l 


dt 


n 


dY 


f    (Y  —     )Y 


+  fo(Y,f3...)Y  -  f(t) 


(37) 


being  the  rules  for  the  procedure  as  follow: 

1.  Obtain  the  Laplace  transform  of  (37)  in  the  usual 
manner  including  initial  conditions. 

2.  Proceed  as  in  the  constant  coefficient  case.   In 

performing  the  long  division  process,  however,  change  the 

value  of  C   •  •  •  C)   at  each  step  in  the  division  process  to 
n      o  ^ 

the  most  recent  values  available  for 
C   =  f  (Y    —     ) 

cn-i  -  fn-i<*.af>---> 


(38) 


Co  ■  V^fp---' 
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In  general 


Ci  =  tl<Y»?jk"'-)  (39) 


In  other  words,  in  the  first  step,  the  values  of  C.  are 

dY 
obtained  from  initial  conditions  for  Y,  tt-.  ...  .   The 

dt 

result  of  the  first  division  is  a  new  value  of  Y  to  be  used 
in  calculating  C  .  C   ,,  ...,  C   for  the  second  step  in  the 
division.   New  values  to  be  used  for  the  derivatives  of  Y 
with  respect  to  t  are  calculated  from  the  new  values  of  the 
dependent  variables.   The  process  is  difficult  to  describe 
in  words,  but  is  simple  in  concept. 

To  illustrate  the  procedure  an  example  is  given: 

Example  3: 


Let        ^   +  Y2  =  1  (40) 

dt 


Carrying  out  the  procedure  outlined  above,  eq.(Ao)  is 
written 


§J  +  C-Y  =  1  (41) 

dt 


The  Laplace  transform,  assuming  C  to  be  constant  and 
expressing  in  powers  of  s    is 


s-2 

Y(s)  =    S   -=  (42) 

1  +  Cs  1 
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Substituting  the  z-forms  of  Table  I  and  dividing  by  T  yields 


Y  *(Z)  =     Tdnoz'1^"2) 

A       (12  +  6CT)  -  24z-1  +  (12-6CT)z" 


Letting  T  =  0.1, 


Y  *(z)  =        0.1  +  l.Oz  -1-  +  z  2 


(43) 


(12  + 0.6C)  -  2  4z  X  +  (12-0.6C)z~2 


Carrying  out  the  long  division  procedure  yields 


00833  +  0.  lOOz  X  +  0.199z-2  +  0.295z"3  +  0.390z~4  +  •  •  •    (44) 
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+  Approximate  solution 


25 


III.   CASE  STUDIES 

A.   INTRODUCTION 

In  order  to  test  the  present  technique,  a  collection 
of  linear  differential  equations  with  varying  coefficients 
were  studied,  whose  basic  format  is: 

Y  +  P(t)Y  =  Q(t) 

where  P(t)  and  Q(t)  are  time  functions.   For  the  purpose  of 
this  work,  nine  equations  in  all  were  processed  successfully 
for  what  P(t)  took  the  values  /t,  t  and  t  ,  and  Q(t)  took 
the  values  l-u(t),  t  and  t  . 

The  collection  of  equations,  the  case  nubmers,  and  the 
sections  in  this  work,  are  as  follows : 


Y  +ftY   =  1 


Case  1 


Y  +   tY  =  1  Case  2      III  B 

Y  +  t2Y  =  1  Case  3 

Y  +ft  Y  =  t  Case  k 

Y  +  tY  =  t  Case  5      HI  C 

Y  +  t2Y  =  t  Case  6 

Y  +  sft  Y  =  t2  Case  7 

Y  +   tY  =  t2  Case  8      III  C 

Y  +  t2Y  =  t2  Case  9 


In  Case  10,  an  example  from  Ref.  2  was  studied,   Comments 
on  this  problem  are  given  in  Case  10  discussion. 
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In  cases  1  to  9  the  same  type  of  procedure  was  followed, 
applying  directly  the  steps  as  in  II  A  of  this  work  and 
experimenting  for  time  intervals  form  0.2  sec.  to  2.0  sec. 
In  order  to  show  the  advantages  of  this  technique  with 
relatively  large  At's,  its  stability  together  with  an 
algorithm  for  computer  solution  and  application  to  hand 
calculation,  i.e.,  small  number  of  division  steps  for  good 
results. 

B.   CASE  l-Y+\/tY=l 

The  equation  Y  +  P(t)Y  =  Q(t)  is  rewritten  for  this 
case  as: 


§+^Y=1       Y(t=o)  ■  ° 


following  the  procedure,  as  in  II  A 


dY 


H  +  PY  =  1  P  =Jt 


the  Laplace  transform  is 
sy  +  Py  =  1/s 


rearranging  and  expressing  in  powers  of  s~  : 


s2y  +  Psy  =  1  y(s2+Ps)  =  1 

.  .    1    .    s-2 


s2  +  Ps    1  +  Ps  -1 
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-k 

Substituting  the   z-form  values   for  s        from  Table   I    (page    15) 


2  -1        -2 

T        1  +  lOz   x  +  z    * 


'(.)    --S2 0.-Z-1)2 


:i+  p  ?1  +  Z"J 


2  1  -  z-1 


Dividing  by  T  and  rearranging: 


J.'W  -   T.lOTz-^Tz-2 ([|5) 

(12  +  6PT)  -  24z  x  +  (12-6PT)z~^ 


where  again  P  is  the  time  varying  coefficient  of  Y  in  the 
differential  equation  and  T  is  the  iteration  time,  or  in 
other  words  the  time  difference  between  results. 

The  following  step  consists  of  dividing  the  denominator 
into  the  numerator  and  substituting  in  P,  its  value  for  the 
corresponding  time  (t  =  nT) .   As  the  exponent  of  z  for  the 
first  value  in  the  quotient  is  zero,  P(t)  =Jt       at  t  =  0 , 
P(t)  =  0  in  the  first  division  step. 

The  value  for  T  first  chosen  was  0.2  sec.  . 

The  algorithm  used  for  computing  the  quotient  (Program  1, 
page  30)  as  the  solution,  for  values  of  T  from  0.2  sec.  up 
to  2  sec.  with  the  "precise"  solution  in  continuous  line  and 
with  the  approximate  solution  as  crosses  follow;  the  "precise" 
solution  was  obtained  by  computer,  using  the  subroutine 
called  INTEG1. 
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Scales  in  the  graphs  (Pigs.  7  to  20)  are  the  same  in 
each  case. 

In  the  numerical  output  tables,  the  time  is  referred  to 
as  T,  the  approximate  solution  is  referred  to  as  Q(N)  and 
the  "precise"  solution  as  PR.   TFIN  is  the  final  time. 

It  can  be  seen  that  the  approximated  solution  even  for 
large  time  iterations  is  very  close  to  the  precise  solution 

FOR  DELTA  T  =  0.20        TFIN  =5.2 

T  Q(N)  PR 

0.0 
0.2 
0.4 
0.6 
0.8 
1.0 
1.2 
1.4 
1.6 
1.8 
2.0 
2.2 
2.4 
2.6 
2.8 
3.0 
3.2 

3.4 
3.6 
3.8 
4.0 
4.2 
4.4 
4.6 
4.8 
5.0 


0.01667 

0.0 

0.19144 

0.19300 

0.36010 

0.36197 

0.49870 

0.50039 

0.60588 

0.60716 

0.68335 

0.68414 

0.73455 

0.73486 

0.76372 

0.76363 

0.77528 

0.77491 

0.77345 

0.77289 

0.76192 

0.76128 

0.74384 

0.74318 

0.72170 

0.72109 

0.69743 

O.69689 

0.67245 

0.67200 

0.64773 

0.64738 

0.62392 

O.62366 

0.60141 

0.60122 

0.58039 

0.58027 

0.56092 

O.56085 

0.54297 

0.54294 

0.52648 

0.52647 

0.51132 

0.51133 

0.49738 

0.49741 

0.48454 

0.48458 

0.47268 

0. 47272 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C       Q(T)=STEP     ,    P(T)=SQRT(T) 
C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

INTEGER    *4ITB(12)/i2*0/ 

REAL    *4RTB(28)/28*0.0/ 

DIMENSION    A(IOO),     Q(IOO),    X(IOO),    XX(IOO),     PR(IOO) 
ITB(3)     =    5 
ITB(4)     =    5 
C 

DO    1    1=1,26 

READ    (5,5)    XX(  I)  ,PR(I ) 

1  CONTINUE 
C 

2  READ    (5,5,END=4)     TD,TFIN 
WRITE     (6,6)     TD,TFIN 
WRITE     (6,7) 

M    =    TFIN/TD 
A(l)    =    TD 
A(2)    =    10.0*TD 
A(3)    =     TD 
A(4)    =    0.0 
C 

DO    3    N=1,M 
T    =    (N-1)*TD 


C 
C 


r 


P    =    SQRT(T) 


Fl    =    12.0+6.0*(P*TD) 

F2    =    -24.0 

F3    =    12.0-6.0*(P*TD) 
C 
C         Q(N)     IS    THE    OUTPUT 


Q(N)     =    A(N)/F1 

A(N+1)     =    A(M+1)-(Q(N)*F2) 

A(N+2)     =    A(N+2)-(Q(N)*F3) 

A(N+3)     =    0.0 

X(N)    =    T 

WRITE     (6,5)     T,Q(N) 

3  CONTINUE 

WRITE    (6,8) 

WRITE    (6,5)     ( XX( I ) ,PR( I) ,1=1,26) 

ITB(l)     =    1 

ITB(2)     =    0 

ITB(12)     =    1 

CALL    DRAWP    ( 26 , X X , PR , I TB, RTB ) 

ITB( 1)     =    3 

ITB(2)     =    2 

CALL    DRAWP    < M , X , Q , I TB, RTB ) 

GO  TO  2 

4  STOP 

5  FORMAT  (2F10.5) 

6  FORMAT  I'IS'FOR  DELTA  T=  •  ,  F6  .  2  ,4X,  ■  TFIN=',F4.1) 

7  FORMAT  (•   T=« , 7X, '  Q(N)=' ) 

8  FORMAT  (•   XX=',7X,'    PR=') 
END 


PROGRAM  1 


30 


x-scale  =1.0  units/inch 
y-scale  =  0.2  units/inch 

FOR  DELTA  T  =  0.2 


FIGURE  7 
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o 
o 


032 


JC3 


x-scale   =1.0  units/inch 


A  T 

-  0.3 

TPIN  =5.2 

T 

Q(N) 

0.0 

0.02500 

0.3 

0.27722 

0.6 

0.49673 

0.9 

0.64695 

1.2 

0.73424 

1.5 

0.77180 

1.8 

0.77421 

2.1 

0.75441 

2.iJ 

0.72252 

2.7 

O.68563 

3.0 

0.64822 

3.3 

0.61282 

3.6 

0.58059 

3.9 

0.55187 

4.2 

0.52653 

4.5 

0.50422 

4.8 

FIGURE  8 

0.48454 
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DELTA   T 

=   0. 

40 

TFIN  =5.2 

T 

Q(N) 

0.0 

0.03333 

0.4 

0.35508 

0.8 

0.60241 

1.2 

0.73386 

1.6 

0.77661 

2.0 

0.76404 

2.4 

0.72370 

2.8 

0.67393 

3.2 

0.62484 

3.6 

0.58087 

4.0 

0.54317 

4.4 

0.5H37 

4.8 

0.48453 

FIGURE 

9 
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FOR  DELTA  T  =  0.1 

>    TFIN  =5.2 

T 

Q(N) 

0.0 

0.4167 

0.5 

0.42489 

1.0 

0.67982 

1.5 

0.77314 

2.0 

0.76570 

2.5 

0.71310 

3.0 

0.64984 

3.5 

0.59170 

4.0 

0.54331 

4.5 

0.50424 

FIGURE    10 
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FOR  DELTA 

T  =  0.6 

TFIN  =5.2 

T 

Q(N) 

0.0 

0.05000 

0.6 

0.48686 

1.2 

0.73288 

1.8 

0.77864 

2.4 

0.72725 

3.0 

0.65099 

3.6 

O.58165 

4.2 

FIGURE  11 

0.52673 
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FOR  DELTA  T  =  0.7 


0.0 
0.7 

1.4 

2.1 

2.8 
4.2 


TFIN  =5.2 

Q(N) 

0.5833 
0.54145 
0.76577 
0.76211 

O.6783I 
0.52680 


FIGURE  12 
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FOR  DELTA 

T  = 

0.8 

TFIN  =5.2 

T 

Q(N) 

0.0 
0.8 
1.6 
2.4 
3.2 
4.0 

0.06667 
0.58920 
0.78249 
0.73260 
0.62874 
0.54381 

FIGURE 

13 
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FOR  DELTA  T  = 


0.0 
0.9 
1.8 

2.7 
3.6 


0.9     TFIN 

Q(N) 

0.07500 

0.63073 
0.78658 
O.69661 
0.58340 


=  5.2 


FIGURE  14 
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FOR  DELTA 

T  = 

■  1 

0 

TFIN  +5.2 

T 

Q(N) 

0.0 
1.0 
2.0 
3.0 
4.0 

0.08333 
O.66667 
0.78105 
0.65849 
0.54411 

FIGURE    15 
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FOR 

DELTA 

T 

0.0 

1.2 
2.4 
3.6 

T 

=  1.2 


TFIN  =5.2 

Q(N) 

0.10000 
0.72408 
0.75053 
0.58590 


FIGURE  16 


40 


n-i'5 


024 


FOR 

DELTA 

T 

0.0 
1.4 
2.8 

T 

=  1.4 


TPIN   =5.2 

Q(N) 

0.11667 
0.76576 
0.70534 


FIGURE   17 
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FOR  DELTA 

T 

=  1. 

6 

TFIN  =5.2 

T 

Q(N) 

0.0 
1.6 

3.2 

0.13333 
0.79526 
0.65424 

FIGURE 

18 
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FOR  DELTA  T 

=1.8     TFIN  =5.2 

T 

Q(N) 

0.0 

1.8 

0.15000 
0.81541 

FIGURE  19 
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DELTA  T  =  2.0 

TFIN  =5.2 

T 

Q(N) 

0.0 
2.0 

0.16667 
0.82843 

FIGURE 

20 
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1.   Case  2  -  Y  +  tY  =  1 

The  calculations  are  the  same  as  In  Case  1  with  the 
only  difference  that  P(t)  is  now  equal  to  t.   So, -when 
dividing  as  in  Equation  (45),  repeated  here  for  convenience, 
the  denominator  into  the  numerator  the  substitution  for  each 
step  is  as  P  =  t . 


Y  *(Z)  =        T  +  lOTz  X  +  Tz~2 

A       (12  +  6PT)  -  24z_1  +  (12-6PT)z~2 


In  the  same  manner,  numerical  and  graphical  output  for 
T  from  0.2  sec.  to  2.0  sec.  (Figs.  21  to  3^4 )  follow. 
Program  2  (page  47 )  shows  a  computer  algorithm  for  solution 


of  this  problem. 


FOR  DELTA  T  =  0.2 


0.0 
0.2 
0.4 
0.6 
0.8 
1.0 
1.2 
1.4 
1.6 
1.8 
2.0 
2.2 
2.4 
2.6 
2.8 
3.0 
3.2 
3.4 
3.6 
3.8 


TFIN  =5.2 

Q(N) 

PR 

O.OI667 

0.0 

0.19608 

0.19735 

0.37707 

0.37933 

0.53018 

0.53293 

0.64664 

0.64931 

0.72264 

0.72478 

0.75927 

0.76058 

0.76154 

0.76198 

0.73700 

0.73670 

0.69413 

0.69335 

0.64099 

0.64000 

0.58425 

0.58329 

0.52880 

0.52801 

0.47769 

0.47713 

0.43241 

0.43208 

0.39333 

0.39319 

0. 36010 

0.36008 

0.33199 

0.33204 

0.30817 

0.30825 

0.28784 

0.28793 
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FOR   DELTA   T 


1  - 

=  0. 

, 2               TFIN 
(Continued) 

=  5. 

.2 

T 

Q(N) 

PR 

4. 

4, 
4. 
4. 
4. 
5. 

.0 
.2 
,4 
,6 
.8 
,0 

0.27033 
0.25507 
0.24162 
0.22966 
0.21893 
0.29023 

0. 
0. 
0, 
0. 
0. 
0. 

,27041 
,25514 
,24168 
,22971 
,21897 
,20926 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C       C(T)=STEP    ,    P(T)=T 

c 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

INTEGER    *4ITB( 12) /12*07 
REAL    *4RTB(28)/28*0.0/ 

DIMENSION    A(IOO),     Q(IOO),    X(IOO),    XX(IOO),     PR(IOO) 
ITB(3)     =    5 
ITB(4)     =    5 
C 

DO    1    1=1,26 

READ    (5,7)    XX(I),PR(I) 

1  CONTINUE 
C 

2  READ    (5,7,END=^)     TD,TFIN 
WRITE    (6,5)     TD,TFIN 
WRITE    (6,6) 

M    =    TFIN/TD 
A(l)    =    TD 
A(2)    =     10.0*TD 
A(3)    =    TD 
A(4)     =    0.0 


C 

c 


r 


DO    3    N=1,M 
T    =    (N-1)*TD 

P    =    T 


Fl    =    12.0+6.0*(P*TD) 

F2    =    -24.0 

F3    =    12.0-6.0*(P*TD) 
C 
C  Q(N)     IS    THE    OUTPUT 


Q(N)     =    A(N)/F1 

A(N+1)     =    A(N+1)-(Q(N)*F2) 

A(Ni-2)     =    A(N+2)-(Q(N)*F3) 

A(N+3)     =    0.0 

X(N)    =     T 

WRITE     (6,7)     T,Q(N) 

3  CONTINUE 

WRITE     (6,8) 

WRITE    (6,7)     (XX( I ) ,PR( I) ,1=1,26) 

ITB(l)     =    1 

ITB(2)     =    0 

ITB(12)     =    1 

CALL    DRAWP    ( 26 , XX , PR , I TB, RTB ) 

ITB(l)     =    3 

ITB(2)     =    2 

CALL  DRAWP  ( M , X , Q , I TB, RTB ) 

GO  TO  2 

4  STOP 

5  FORMAT  CI', 'FOR  DELTA  T= ' , F6. 2 , 4X, «  TFIN=',F4.1) 

6  FORMAT  (•   T=,,7X,«  Q(N)=') 

7  FORMAT  (2F10.5) 

8  FORMAT  (•   XX=' ,7X,'    PR=« ) 
END 


PROGRAM  2 


^7 


x-scale  =  1.0  units/inch 
y-scale  =0.2  units/inch 

FOR  DELTA  T  =  0.2 


FIGURE  21 
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FOR  DELTA  T 

=  0.3 

TFIN  =5.2 

T 

QCN) 

0.0 

0.02500 

0.3 

0.28708 

0.6 

0.52675 

0.9 

0.68665 

1.2 

0.75759 

1.5 

0.75202 

1.8 

0.69513 

2.1 

0.61402 

2.H 

0.52986 

2.7 

0. 45^88 

3.0 

0.39355 

3.3 

0.3^5 

3.6 

0.30809 

3.9 

O.27869 

4.2 

0.25500 

4.5 

0.235^3 

4.8 

0.21891 

FIGURE 
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uu 
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FOR  DELTA 

T   = 

0. 

40 

TFIN   = 

T 

Q(N) 

0.0 

0.03333 

0.4 

0.37037 

0.8 

0.63857 

1.2 

0.75516 

1.6 

0.73782 

2.0 

0.64408 

2.4 

0.53138 

2.8 

0.43354 

3.2 

0. 36022 

3.6 

0.30795 

4.0 

0.27013 

4.4 

0.24150 

4.8 

0.21887 

FIGURE    2  3 
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FOR  DELTA 

T   =    0 

5 

TFIN   =5.2 

T 

Q(N) 

0.0 
0.5 

0.04167 
0.44444 

1.0 

0.71H1 

1.5 
2.0 

0.75151 
0.64646 

2.5 
3.0 

0.50660 
0.39427 

3-5 
1.0 

0.31923 
0.26995 

4.5 

0.23529 

FIGURE 

24 
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FOR  DELTA  T 

=  0.6 

TFIN  =5.2 

T 

Q(N) 

0.0 

0.5000 

0.6 

0.50847 

1.2 

0.74776 

1.8 

0.70037 

2.4 

0.53615 

3.0 

0.39480 

3-6 

0.30744 

4.2 

0.25460 

FIGURE 

25 
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FOR  DELTA 

T  = 

0.7 

TFIN  =5.2 

T 

Q(N) 

0.0 

0.05833 

0.7 

0.56225 

1.4 

0.75470 

2.1 

0.62530 

2.8 

0.43722 

3.5 

0.31854 

4.2 

0.25438 

FIGURE 

26 
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FOR  DELTA  T 

T 

0.0 
0.8 
1.6 
2.4 
3.2 
4.0 


=  0.8 


TPIN  =5.2 

Q(N) 

0.06667 
0.60606 
0.73910 
0.54392 
0.36042 
0.26888 


FIGURE  2  7 


54 


1-. 

002 

JC3                          0Q? 

FOR   DELTA 

T 

=    0 

9 

TFIN   =5.2 

T 

Q(N) 

0.0 
0.9 
1.8 
2.7 
3.6 

0.07500 
0.64057 
0.70781 
0.46704 
0.30519 

FIGURE 

28 
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FOR  DELTA  T 

=  1 

0 

TFIN  =5.2 

T 

Q(N) 

0.0 
1.0 
2.0 
3-0 

0.08333 
0.66667 
0.66667 
0.40000 
0.26667 

FIGURE 

29 

56 


FOR  DELTA  T 

T 

0.0 
1.2 
2.4 
3.6 


=  1.2 


TPIN  =5.2 

Q(N) 
0.10000 

0.69767 
0.57186 
0.30012 


FIGURE  30 
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FOR  DELTA  T 

T 

0.0 
l.k 

2.8 


=  l.H 


TFIN  =5.2 

Q(N) 

0.11667 
0.70707 
0.47775 


FIGURE  31 
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FOR  DELTA  T 

T 

0.0 
1.6 
3.2 


=1.6    TFIN  =5.2 

Q(N) 

0.13333 
0.70175 
0.39^24 


FIGURE  32 
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ELTA 

T 

=  1.8 

TFIN  =  5 

2 

T 

Q(N) 

0.0 
1.8 

FIGURE 

33 

0.15000 
0.68702 

60 


FOR  DELTA  T 

T 

0.0 
2.0 


=  2.0 


TFIN  =5.2 

Q(N) 

0.16667 
O.66667 


FIGURE  34 
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2.   Case  3  -  Y  +  t2Y  =  1 

As  in  Case  1  and  Case  2,   the  same  calculations  are 
followed  and  only  in  the  division,  the  substitution  of 
P(t)  =  t2  is  the  sole  difference. 

Output  for  several  time  iterations  follows  as  in  the 
previous  cases  (Figs.  35  to  48).   Computer  algorithm  is  shown 
in  Program  3  (page  63). 

FOR  DELTA  T  =  0.2  TFIN  =5.2 

T  Q(N)         PR 

0.0 
0.2 
0.4 
0.6 
0.8 
1.0 
1.2 
1.4 
1.6 
1.8 
2.0 
2.2 
2.4 
2.6 
2.8 
3.0 
3.2 
3.4 
3.6 
3.8 
4.0 
4.2 
4.4 
4.6 
4.8 
5.0 

Equation  (45)  Is  repeated  here  for  convenience: 


~,     ,  T  +  10Tz  -1  +  Tz~2 


0.01667 

0.0 

0.19920 

0.19960 

0.39213 

0.39366 

0.56550 

0.56858 

0.70032 

0.70472 

0.77773 

0.78240 

0.78667 

0.79011 

0.73026 

0.73132 

0.62669 

0.62542 

0.50321 

0.50080 

0.38584 

0.38376 

0.29077 

0.28980 

0.22210 

0.22208 

0.17552 

0.17589 

0.14398 

0.14433 

0.12163 

0.12183 

0.10482 

0.10492 

0.09160 

0.09165 

0.08088 

0.08091 

0.07204 

0.07205 

0.06462 

0.06464 

0.05833 

0.05834 

0.05294 

0.05295 

0.04828 

0.04830 

0.04422 

0.04423 

0.04067 

0.04067 

(12  +  6PT)  -  24z_1  +  (l2-6PT)z" 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C       Q(T)=STEP     ,    P(T)=T**2 
C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

INTEGER    *4ITB<12)/12*0/ 

REAL    *4RTB(28)/28*0.0/ 

DIMENSION    A(IOO),     Q(IOO),    X(IOO),    XX(IOO),     PR(IOO) 
ITB(3)     =    5 
ITB(4)     =    5 
C 

DO    1    1=1,26 

READ    (5,7)    XX(  I)  ,PRU  ) 

1  CONTINUE 
C 

2  READ    (5,7,END=4)     TD,TFIN 
WRITE     (6,5)     TD,TFIN 
WRITE     (6,6) 

M    =    TFIN/TD 
A(l)    =    TD 
A(2)    =     10.0*TD 
A(3)    =    TD 
A(4)     =    0.0 
C 

DO    3    N=1,M 
T    =    (N-1)*TD 


C 
C 


c 


P    =    T**2 


Fl  =    12.0+6.0*(P*TD) 

F2  =    -24.0 

F3  =    12.0-6.0*(P*TD) 
C 

C  Q(N)  IS    THE    OUTPUT 


Q(N)     =    A(N)/F1 

A(N+1)     =    A(N+1)-(Q(N)*F2) 

A(N+2)     =    A(N+2)-(Q(N)*F3) 

A(N+3)     =    0.0 

X(N)    =    T 

WRITE     (6,7)    T,Q(N) 

3  CONTINUE 

WRITE    (6,8) 

WRITE     (6,7)     (XX( I ) , PR( I) , 1=1,26) 

ITB(l)     =    1 

ITB(2)     =    0 

ITB(12)     =    1 

CALL    DRAWP    ( 26 , XX , PR ,  ITB, RT B ) 

ITB( 1)     =    3 

ITB(2)     =    2 

CALL    DRAWP    ( M , X , Q , I TB, RTB ) 

GO    TO    2 

4  STOP 

5  FORMAT  CI', 'FOR  DELTA  T= «  ,  F6 . 2 , 4X,  •  TFIN=«,F4.1) 

6  FORMAT  ( '   T  =  '  ,7X,  •  Q(N)='  ) 

7  FORMAT  (2F10.5) 

8  FORMAT  ('   XX=,,7X,'    PR=») 

END 


PROGRAM  3 
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x-scale  =1.0  units/inch 
y-scale  =  0.2  units/inch 


FOR  DELTA  T  =  0.2 


FIGURE  35 


6'l 


11 
a 
o 


CI 

a 


a  L 


,-»1  1 

L.U-1 

0G2 

)C3                           03  *, 

FOR  DELTA 

T    = 

0. 

3 

TFIN   =5.2 

T 

Q(N) 

0.0 

0.02500 

0.3 

0.29600 

0.6 

0.56168 

0.9 

0.74128 

1.2 

0.78225 

1.5 

0.68283 

1.8 

0.50631 

2.1 

0.33719 

2.1 

0.22218 

2.7 

0.15773 

3.0 

0.12138 

3.3 

0.09778 

3.6 

0.08086 

3.9 

O.O6816 

it. 2 

O.05832 

1.5 

0.05053 

1.8 

0.01122 

FIGURE 

36 

65 


FOR  DELTA 

T  =  0.4 

TFIN  =5.2 

T 

Q(N) 

0.0 

0.03333 

0.4 

0.38760 

0.8 

0.68723 

1.2 

0.77582 

1.6 

O.62989 

2.0 

0.39299 

2.4 

0.22240 

2.8 

0.14260 

3.2 

0.10466 

3.6 

0.08082 

ii.Q 

0.06460 

4.4 

0.05293 

4.8 

0.04422 

FIGURE 

37 
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FOR  DELTA  T  = 

0.5 

TFIN  =5.2 

T 

Q(N) 

0.0 

0.04167 

0.5 

0.47059 

1.0 

0.75294 

1.5 

0.68141 

2.0 

0.39906 

2.5 

0.19512 

3.0 

0.12007 

3.5 

O.O8613 

4.0 

0.06447 

4.5 

0.05057 

FIGURE 

38 
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DELTA 

T  =  0.6 

TFIN  =5.2 

T 

Q(N) 

0.0 

0.05 

0.6 

0.54152 

1.2 

0.75631 

1.8 

• 

0.52210 

2.4 

0.22530 

3.0 

0.11783 

3.6 

0.08177 

4.2 

0.05783 

FIGURE 

39 
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FOR  DELTA  T 

=  0 

7 

TFIN  =5.2 

T 

Q(N) 

0.0 

0.05833 

0.7 

0.59752 

1.4 

0.70881 

2.1 

0.36272 

2.8 

0.13431 

3.5 

0.08809 

4.2 

0;05721 

FIGURE 

40 
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U1 

002 

003                           00: 

FOR  DELTA 

T 

=    0 

8 

TFIN  =5.2 

T 

Q(N) 

0.0 
0.8 
1.6 
2.4 
3.2 
4.0 

0.06667 
0.63694 
0.62939 
0.23756 
0.09620 
O.O6786 

FIGURE 

41 
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FOR  DELTA  T 

T 

0.0 

0.9 
1.8 

2.7 
3.6 


0.9    TPIN  =5.2 

Q(N) 

0.07500 
0.65958 
0.53668 
0.15283 
0.08072 


FIGURE  iJ2 
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FOR  DELTA   T 

T 

0.0 
1.0 
2.0 
3.0 
4.0 


=    1.0 


TFIN   =5.2 

Q(N) 

0.08333 
O.66667 
0. 114444 

0.10101 
O.07183 


FIGURE  ^3 
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FOR  DELTA 

T 

=  1 

.2 

TFIN  =5.2 

T 

Q(N) 

0.0 
1.2 
2.4 
3.6 

0.10000 
0.64378 
0.28895 
0.05587 

FIGURE 

44 
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"1  -1 

CGI                          0G2 

FOR  DELTA   T   =    1.4 

TFIN   =5.2 

T 

Q(N) 

0.0 
1.4 
2.8 

0.11667 
0.59022 
0.18194 

FIGURE 

45 
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FOR  DELTA  T 

T 

0.0 
1.6 
3.2 


1.6    TFIN  =5.2 

Q(N) 

0.13333 
0.52493 
0.11422 


FIGURE  46 
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FOR  DELTA  T 

T 

0.0 
1.8 


=  1.8 


TPIN  =5.2 

Q(N) 

0.15000 
0.45965 


FIGURE  47 
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FOR  DELTA 

T 

=  2 

.0 

TFIN  =5.2 

T 

Q(N) 

0.0 
2.0 

0.16667 
0.40000 

FIGURE 

48 
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C.   CASE  4-Y+V^tY  =  t 

For  this  case,  in  the  equation  Y  +  P(t)Y  =  Q(t), 
P(t)  =/t   and  Q(t)  =  t,  so  following  the  procedure  the 
equation  is  rewritten  as: 


f  +  PY  -  *        Y(t=o)  -  ° 


where   P  =  J~t      .   The  Laplace  transform  is 

.2 


so 


sy  +  Py  =  l/s' 


s3y  +  Ps2y  =  1  y(s3  +  Ps2)  =  1 


y  = 


1        s-3 


s3  +  Ps2    1  +  Ps  l 


Substituting  the  z-forms  of  s~   from  Table  I  one  gets 


T3  z_1 +  z~2 


..(.,  .   2  (i-.-l)3 


T  1  +  z-1 
1  +  P  £       , 

2  1  -  z"1 


Dividing  by  T  and  ordering  in  crescent  powers  of  z~ 
in  numerator  and  denominator: 


2-1     2  -? 
*r    ^  -            T  Z    +  T  z  c  fU(-s 

yA*(z)  =  ; rp ir   (*»6) 

R  (2+Pt)  -  (6+PT)z  i  +  (6-PT)z  ^  +  (PT-2)z  j 
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Now,  substituting  for  the  time  iteration  wanted  and 
replacing  P  =  vt   for  its  value  at  the  division  step,  the 
solution  is  obtained. 

In  the  same  way  as  before,  plots  and  outputs  follow 
for  time  iterations  from  0.2  sec.  to  2.0  sec.  (Figs.  49  to 
60). 

Algorithm  for  computer  solution  of  this  problem  is  shown 
in  Program  4  (page  8o) • 


FOR  DELTA  T  =  0.2 

TFIN  =5.2 

T 

Q(N) 

PR 

0.2 

0.01914 

0.0 

0.4 

0.07363 

0.01949 

0.6 

0.15683 

0.07449 

0.8 

0.26131 

0.15817 

1.0 

0.37994 

0.26299 

1.2 

0.50647 

0.38181 

1.4 

0.63576 

0.50837 

1.6 

0.76391 

0.63756 

1.8 

O.88813 

0.76551 

2.0 

1.00662 

0.88948 

2.2 

I.II838 

1.00770 

2.4 

1.22302 

1.11920 

2.6 

1.32060 

1.22360 

2.8 

1.41147 

1.32100 

3.0 

1.49613 

1.41170 

3.2 

1.57520 

1.49620 

3.4 

1.64928 

1.57520 

3.6 

1.71898 

1.64930 

3.8 

1.78485 

1.71890 

4.0 

1.84739 

1.78480 

4.2 

1.90703 

1.84730 

4.4 

1.96413 

1.90700 

4.6 

2.01902 

1.96410 

4.8 

2.07196 

2.01900 

5.0 

2.12316 

2.07200 

5.2 

2.17281 

2.12320 
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CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC3C 

c 
c 
c 
ccc 


C(T)=RAMP       £       P(T)=SQRT(T) 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

INTEGER  *4ITB<12)/12*0/ 

REAL  *4RTB<28)/28*0.0/ 

DIMENSION    A(100),     Q( 1 00} ,    X ( 100 ) ,    XX(IOO),     PR(IOO) 

ITB(3)     =    5 

ITB(*)     =    5 

WRITE    16,8) 

DO    1    1=1,26 

READ    (5,7)    XX( I)  , PR( I ) 

WRITE     (6,7)    XX(I),PR(I) 

1  CONTINUE 

2  READ    (5,7,END=<+)     TD.TFIN 
WRITE    (6,5)     TD,TFIN 
WRITE     (6,6) 

X(l)    =    0.0 
M    =    TFIN/TD 
A(l)    =    TD**2 
A(2)    =    TD**2 
A(3)    =    0.0 
A(4)    =    0.0 


C 
C 


DO    3    N=1,M 
T    =    N*TD 
P    =    SQRT(T) 
Fl    =    2.0+P*TD 
F2    =    -( 6.0+P*TD) 
F3    =    6. 0-P*TD 
F4   =    (P*TD)-2.0 

Q(N)    =    A(N)/F1 


A(N+1) 

A(N+2) 

A(N+3) 

A(N+4) 

WRITE 

X(N)    = 


A(N+1)-(Q(N)*F2) 

A(N+2)-( Q(N)*F3) 

A(N+3)-(Q(N)*F4) 

0.0 

7)     T,Q(N) 


3    CONTINUE 


ITB(l)     =    1 

ITB(2)     =    0 

RTB(6)     =    TD 

ITB(12)     =    1 

CALL    DRAWP    ( 26, XX , PR, I TB, RTB ) 

ITB(l)     =    3 

ITB(2)     =    2 

CALL  DRAWP  ( M , X , Q , I TB, RTB ) 

GO  TO  2 

STOP 


FORMAT 
FORMAT 

FORMAT 
FORMAT 
END 


CI' 
(• 

(2F10.5) 
(•   XX=« 


FOR  DELTA  T=« ,F6, 
Q(N)=' ) 

7X,'    PR='  ) 


2,4X, •TFIN=' ,F4.1 ) 


PROGRAM  k 
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+ 


x-scale  =1.0  units/inch 
y-scale  =  0.5  units/inch 

FOR  DELTA  T  =  0.2 


FIGURE  49 
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qz: 

CZ2 

«J»-W 

FOR  DELTA 

T  =    0.3 

TFIN   =5.2 

T 

Q(N) 

0.3 

0.04158 

0.6 

0.15514 

0.9 

0.31700 

1.2 

0.50407 

1.5 

0.69801 

1.8 

0.88639 

2.1 

1.06213' 

2.4 

1.22223 

2.7 

1.36641 

3.0 

1.49595 

3.3 

1.61280 

3.6 

1.71905 

3.9 

1.81664 

4.2 

1.90718 

4.5 

1.99200 

4.8 

2.07212 

5.1 

2.14834 

FIGURE 

50 
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4- 


FOR  DELTA 

T  =  0.4 

TFIN  =5.2 

T 

Q(N) 

0.4 

0.07102 

0.8 

0.25620 

1.2 

0.50068 

1.6 

0.75898 

2.0 

1.00322 

2.4 

1.22111 

2.8 

1.41066 

3.2 

1.57505 

3.6 

1.71916 

4.0 

1.84769 

4.4 

1.96446 

4.8 

2.07227 

5.2 

FIGURE  51 

2.17312 
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DELTA 

T  = 

0.5 

TFIN  =5-2 

T 

Q(N) 

0.5 

0.10622 

1.0 

0.36996 

1.5 

0.69092 

2.0 

1.00060 

2.5 

1.26987 

3.0 

1.49539 

3.5 

1.68484 

H.O 

1.84788 

4.5 

1.99233 

5.0 

2.12360 

FIGURE 

52 
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->Z 

f]  "•  1 

CQ2 

a 

FOR   DELTA 

T   =    0 

6 

TFIN   =5.2 

T 

Q(N) 

0.6 
1.2 

1.8 
2.4 
3.0 
3.6 
4.2 
4.8 

0.14606 
0.49082 

O.87667 
1.21782 
1.49500 

1.71944 

1.90777 
2.07258 

FIGURE 

53 
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FOR  DELTA 

T  = 

0. 

7 

TFIN  =5.2 

T 

Q(N) 

0.7 

0.18951 

1.4 

0.61452 

2.1 

1.05164 

2.8 

1.40840 

3.5 

1.68514 

4.2 

1.90808 

4.9 

2.09855 

FIGURE 

54 
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FOR  DELTA  T  = 


0.8 
1.6 

2.4 
3.2 
4.0 
4.8 


0.8     TFIN  =5.2 

Q(N) 

0.23568 
0.73797 
1.21295 
1.57461 
1.84884 
2.07304 


FIGURE  55 
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FOR  DELTA  T  =  0.9 


TFIN 


5.2 


T 


0.9 
1.8 

2.7 
3.6 

4.5 


Q(N) 

0.28383 
0.85903 
1.35987 
1.72042 
1.99350 


FIGURE  56 
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FOR  DELTA  T 

T 

1.0 

2.0 
3.0 
4.0 

5.0 


=1.0    TFIN  =5.2 

Q(N) 

0.33333 
0.97631 
1.49299 
1.85001 
2.12461 


FIGURE  57 
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FOR  DELTA 

T  =  1.2 

TFIN  =5.2 

T 

Q(N) 

1.2 

2.4 
3.6 
4.8 

0.43445 
1.19662 

1.72293 
2.07450 

FIGURE 

58 

90 


FOR  DELTA 

T 

=  1.4 

TFIN  =5.2 

T 

Q(N) 

1.4 
2.8 
4.2 

0.53603 
1.39641 
1.91440 

FIGURE 

59 
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FOR  DELTA  T 

T 

1.8 
3.6 


=  1.8 


TFIN  =5.2 

Q(N) 

0.73387 
1.73869 


FIGURE  60 
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FOR  DELTA  T  =  2. 

T 
2.0 


TFIN  =5.2 

Q(N) 
0.82843 


FIGURE  61 
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1.   Case  5-Y+tY  =  t 

As  case  2  follows  case  1,  as  calculations  are 
concerned,  case  5  follows  case  4  with  the  only  difference 
that  P(t)  is  now  equal  to  t  and  this  must  be  the  value 
used  to  substitute  in  P  before  and  during  the  division  steps, 

Graphs,  plots  and  solution  values  for  several  time 
iterations  follow  (Pigs.  62-75).   Computer  algorithm  for 
solution  of  this  problem  is  shown  in  Program  5  (page  95  ) . 


FOR  DELTA  T  =  0.2 


T 


0.0 
0.2 
0.4 
0.6 
0.8 
1.0 
1.2 
1.4 
1.6 
1.8 
2.0 


2, 
2. 
2. 
2, 
3. 
3. 


TFIN  =5.2 

Q(N) 

PR 

0.0 

0.0 

0.01961 

0.01980 

0.07617 

0.07688 

0.16332 

0.16473 

0.27178 

0.27385 

0.39094 

0.39347 

0.51058 

0.51324 

0.62220 

0.62468 

0.71991 

0.72186 

0.80061 

0.80209 

0.86375 

0.86466 

0.91066 

0.91107 

0.94380 

0.94386 

0.96617 

0.96595 

0.98040 

O.98016 

0.98914 

O.98889 

0.99424 

0.99402 

0.99707 

0.99691 

0.99857 

0.99847 

0.99933 

0.99927 

0.99970 

0.99966 

0.99986 

0.99985 

0.99993 

0.99994 

0.99996 

0.99997 

0.99997 

0.99999 

0.99997 

1.00000 

3.4 
3.6 
3.8 
4.0 
4.2 
4.4 
4.6 
4.8 
5.0 

Equation  (46)  is  repeated  here  for  convenience 
Y  »(z)  = 


m2  -1  ,  ^2  -2 
T  z    +  T  z 


(2+PT)  -  (6+PT)z-1  +  (6-PT)z  2  +  (PT-2)z  3 
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CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC3C 

c 

C       QIT)=RAMP       &       P(T)=T 
C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

INTEGER    *4ITB( 12)/12*0/ 
REAL    *4RTB(28)/28*0.0/ 

DIMENSION    A(IOO),     Q(IOO),     X(IOO),    XX(IOO),     PR(IOO) 
ITB(3)     =    5 
ITB(4)     =    5 
WRITE    (6,8) 
C 

DO    1    1=1,26 

READ    (5,7)    XX(I) ,PR(I ) 

WRITE     (6,7)     XX(I),PR(I) 

1  CONTINUE 
C 

2  READ    (5,7,END=4)     TD,TFIN 
WRITE    (6,5)     TD,TFIN 
WRITE    (6,6) 

X(l)    =    0.0 
M    =    TFIN/TD 
A(l)    =    TD**2 
A(2)    =    TD**2 
A(3)    =    0.0 
A<4)    =    0.0 
C 

DO    3    N=1,M 

T    =    N*TD 

P    =    T 

Fl    =    2.0+P*TD 

F2   =    -( 6.0+P*T0) 

F3    =    6.0-P*T0 

F4    =    (P*TD)-2.0 

r 

Q(N)    =    A(N)/F1 
C 

A(N+1)     =    A(N+1)-(Q(N)*F2) 
A(N+2)     =    A(N+2)-(Q(N)*F3) 
A(N+3)     =    A(N+3)-(Q(N)*F4) 
A(N+4)     =0.0 
WRITE     (6,7)     T,Q(N) 
X(N)    =    T 

3  CONTINUE 
C 

ITB(l)     =    1 

ITB(2)     =    0 

RTB(6)     =    TD 

ITB(12)     =    1 

CALL    DRAWP    ( 26 , XX , PR , I TB, RT B ) 

ITB(l)     =    3 

ITB(2)     =    2 

CALL    DRAWP    (  M  ,X  ,  Q,  I  TB,  RTB  ) 

GO    TO    2 

4  STOP 
C 

5  FORMAT  ('IS'  FOR  DELTA  T= ' , F6. 2 , 4X , ' TF I N=' ,F4. 1 ) 

6  FORMAT  ( •   T  =  «  ,7X,  •  Q(N)  =  »  ) 

7  FORMAT  (2F10.5) 

8  FORMAT  (•   XX=',7X,«    PR=«) 
END 


PROGRAM  5 
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I  I  I  I  I  I  I-  + 


x-scale 
y-scale 


1.0  units/inch 
0.2  units/inch 


DELTA  T  =  0.2 


FIGURE  62 
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o 


H H 1 1 1 1 —  + 


002 


00  i 


)ELTA  T  =  0 .  3 

TPIN  =5.2 

T 

Q(N) 

0.3 

0.04306 

0.6 

0.16158 

0.9 

0.32779 

1.2 

0.50723 

1.5 

0.67015 

1.8 

0.79871 

2.1 

0.88826 

2.1 

0.9^372 

2.7 

0.97436 

3.0 

0.98948 

3.3 

0.99613 

3.6 

0.99873 

3.9 

0.99963 

1.2 

0.99990 

1.5 

0.99997 

1.8 

0.99999 

5.1 

0.99999 

FIGURE  63 
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+■ 


032 

0C3 

FOR  DELTA 

T    = 

0.4 

TFIN   =5.2 

T 

Q(N) 

0.4 

0.07407 

0.8 

0.26564 

1.2 

0.50253 

1.6 

0.71358 

2.0 

0.86088 

2.4 

0.94360 

2.8 

0.98120 

3.2 

0.99496 

3.6 

0.99895 

4.0 

0.99984 

4.4 

0.99998 

4.8 

1.00000 

5.2 

1.00000 

FIGURE    6  4 
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DELTA 

T  =  0 

5 

TFIN  = 

T 

Q(N) 

0.5 

0.11111 

1.0 

0.37776 

1.5 

0.66061 

2.0 

0.85859 

2.5 

0.95649 

3.0 

0.99067 

3.5 

0.99876 

4.0 

0.99992 

4.5 

1.00000 

5.0 

1.00000 

FIGURE 

65 

5.2 
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FOR  DELTA  T  =  0.6     TFIN  =  5.2 


6 

2 
8 
4 


3.0 
3.6 

4.8 


Q(N) 

0.1525^ 
0.48903 
0.78765 
0.94321 
0.99163 
0.99960 
1.00001 
1.00000 


FIGURE  66 


100 


DELTA 

T  =  0 

7 

TFIN  =5.2 

T 

Q(N) 

0.7 

0.19679 

l.ij 

0.59300 

2.1 

O.88036 

2.8 

0.98399 

3.5 

0.99986 

4.2 

1.00001 

4.9 

1.00000 

FIGURE 

67 

101 


DELTA  T  =  0 . 8 

TFIN  =5.2 

T 

Q(N) 

0.8 
1.6 
2.4 
3.2 
4.0 
4.8 

0.24242 
0.68588 
0.94231 
0.99899 
1.00011 
0.99998 

FIGURE 

68 

102 


FOR  DELTA 

T  =  0 

9 

TFIN  =5.2 

T 

Q(N) 

0.9 

1.8 

2.7 

3.6 

4.5 

0.28826 
0.76603 
0.97993 
1.00165 
0.99966 

FIGURE 

69 

103 


FOR  DELTA 

T  = 

1. 

0 

TFIN  =5.2 

T 

Q(N) 

1.0 
2.0 
3.0 
4.0 
5.0 

0.33333 
0.83333 
1.00000 
1.00000 
1.00000 

FIGURE 

70 
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FOR  DELTA 

T 

=  1 

2 

TFIN  =5.2 

T 

Q(N) 

1.2 
2.4 
3.6 
4.8 

0.41860 
0.93328 
1.00929 
0.99722 

FIGURE 

71 

105 


FOR  DELTA  T 

T 

1.4 

2.8 
4.2 


=  1.4 


TFIN  =5.2 
Q(N) 

0.49495 
0.99659 
1.00083 


FIGURE  72 


106 


FOR  DELTA 

T  = 

l.( 

TFIN  =5.2 

T 

Q(N) 

1.6 
3.2 
4.8 

0.56140 
1.03450 
0.98888 

FIGURE 

73 
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+ 


FOR  DELTA  T  =  1.8 

TFIN  =5.2 

T 

Q(N) 

1.8 
3.6 

0.61832 
1.05581 

FIGURE  7*J 

108 


FOR  DELTA  T  =  2.0 

TFIN  =5.2 

T 

Q(N) 

2.0 

0.66667 
1.06667 

FIGURE  75 
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2.   Case  6  -  Y  +  t  Y  =  t 

Calculations  follow  exactly  Case  4  except  for  ■ 

2 
P(t)  =  t  .   In  the  division  steps,  substitution  of  P  for 

2 
Its  value  t  must  be  done  in  this  case. 

Solutions  for  several  time  iterations  follow 

(Figs.  76  to  89).   Computer  solution  is  in  Program  6  (page 

111). 


FOR  DELTA  T  =  0.2 


TFIN  -  5.2 


Q(N) 


0.2 
0.4 
0.6 
0.8 
1.0 
1.2 
1.4 
1.6 
1.8 
2.0 


2 
2 

2 
2 


3.0 
3.2 
3.4 
3.6 
3.8 
4.0 


4 
4 
4, 
4 


5.0 
5.2 


01992 
07858 
17116 
28666 
40755 
0.51294 
0.58451 
0.61301 
0.60127 
0.56176 
51015 
45891 
41443 
37796 
34824 
32354 
30253 
0.28433 
0.26835 
0.25418 
0.24150 
0.23008 
0.21973 
0.21029 
0.20166 
0.19372 


PR 


0.01997 
0.07898 
0.17243 

0.28923 
0.41144 
0.51742 
0.58837 
0.61523 
0.60163 
56095 
50914 

45829 
41425 
37801 
34833 
32360 
30257 
0.28436 

0.26837 
0.25420 
0.24152 
0.23009 
0.21974 
0.21030 
0.20167 


A 


Equation  (46)  is  repeated  here  for  convenience 
»(z) 


?  -]     2-2 

T^z  L    +  T  z 


(2  +  PT)  -  (6  +  PT)z  ]-+(6-PT)z  2  +  (PT-2)z  3 
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CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC3C 

c 

C       C(T)=RAMP       £       P(T)=T**2 
C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

INTEGER    *4ITB(12)/12*0/ 
REAL    *4RTB(28)/28*0.0/ 

DIMENSION    A(IOO),     Q(IOO),    X(IOO),    XX(IOO),     PR(IOO) 
ITB(3)     =    5 
ITB(4)     =    5 
WRITE     (6»8) 
C 

DO    1    1=1,26 

READ    (5,7)     XX( I) ,PR(I) 

WRITE     (6,7)     XX(  I)  ,PR(  I ) 

1  CONTINUE 
C 

2  READ    (5,7,END=^)     TD,TFIN 
WRITE    (6,5)     TD,TFIN 
WRITE     (6,6) 

X(l)    =    0.0 

M    =    TFIN/TD 

A(l)    =    TD**2 

A(2)    =    TD**2 

A(3)    =    0.0 

A(4)    =    0.0 
C 

DO    3    N=1,M 

T    =    N*TD 

P    =    T**2 

Fl    =    2.0+P*TD 

F2    =    -( 6.0+P*TD) 

F3    =    6.  0-P*TD 

F4    =    (P*TD)-2.0 
C 

Q(N)    =    A(N)/F1 
C 

A(N+1)     =    A(N+1)-(Q(N)*F2) 

A(N+2)     =    A(N+2)-(Q(N)*F3) 

A(N+3)     =    A(N+3)-(Q(N)*F4) 

A(N+4)     =    0.0 

WRITE    (6,7)    T,&(N) 

X(N)    =    T 

3  CONTINUE 
C 

ITB(l)     =    1 

ITB(2)     =    0 

RTB(6)     =    TD 

ITB(12)     =    1 

CALL    DRAWP    ( 26 , XX , PR  ,  I TB, RTB ) 

ITBd  )     =    3 

ITB(2)     =    2 

CALL    DRAWP    (  M  ,  X  ,  Q  ,  I  TB,  R  TB  ) 

GO    TO    2 

4  STOP 
C 

5  FORMAT     CI',1     FOR    DELTA    T  =  '  ,  F6.  2  ,4X  ,  •  TF  I  N=  •  ,F4. 1  ) 

6  FORMAT     ( •       !  =  •  .  7X,  '     Q(N)  =  '  ) 

7  FORMAT     (2F10.5) 

8  FORMAT     ('       XX=,,7X,«  PR=») 
END 


PROGRAM    6 
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+ 


x-scale  =1.0  units/Inch 
y-scale  =0.2  units/inch 

FOR  DELTA  T  =  0.2 

FIGURE  76 
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FOR  DELTA  T  =  0 .  3 


TPIN  =5.2 


Q(N) 


04440 
16964 
34372 
50736 
60020 
0.60070 
0.53792 
45981 
39529 
34810 

31257 
28430 
0.26104 
0.24149 
0.22477 
0.21028 
0.19760 


FIGURE  77 
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+ 


DELTA 

T   = 

OJ 

\          TFIN   = 

T 

Q(N) 

0.4 

0.07752 

0.8 

0.27929 

1.2 

0.49964 

1.6 

0.60565 

2.0 

0.56420 

2.4 

0.46136 

2.8 

0.37768 

3.2 

0.32332 

3.6 

0.28429 

4.0 

0.25415 

4.4 

0.23007 

4.8 

0.21028 

5.2 

0.19371 

FIGURE 

78 

5.2 


11' 


DELTA  T  =  0.5 

TFIN  =5.2 

T 

Q(N) 

0.5 

0.11765 

1.0 

0.38824 

1.5 

0.58635 

2.0 

0.56576 

3.0 

0.34709 

3.5 

0.29320 

4.0 

0.25405 

1.5 

0.22480 

5.0 

0.20162 

FIGURE 

79 

r 

Q2 

DELTA 

T    = 

0.6 

TFIN    =5.2 

T 

Q(N) 

0.6 

0.16245 

1.2 

0.47829 

1.8 

0.59415 

2.4 

0.46797 

3.0 

0.34576 

3.6 

0.28482 

4.2 

0.24117 

4.8 

0.21043 

FIGURE    80 
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FOR  DELTA 

T  = 

0 

7 

TFIN  =5.2 

T 

Q(N) 

0.7 

0.20913 

1.4 

0.53871 

2.1 

0.54812 

2.8 

0.37850 

3.5 

0.29218 

4.2 

0.24177 

4.9 

0.20568 

FIGURE 

81 
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FOR  DELTA 

T    = 

OJ 

3 

TFIN   =5.2 

T 

Q(N) 

0.8 
1.6 

2.4 
3.2 
4.0 
4.8 

0.25478 
0.56796 
0.48014 
0.31670 
O.25669 
0.20888 

FIGURE 

82 
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FOR  DELTA  T 

T 

0.9 
1.8 
2.7 

3.6 

4.5 


0.9    TFIN  =5.2 

Q(N) 

0.29681 
0.57104 
0.41198 
0.27744 
0.22788 


FIGURE  83 
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FOR  DELTA  T 

T 

1.0 
2.0 
3.0 
4.0 
5.0 


=  1.0 


TFIN  =5.2 

Q(N) 

0.33333 
0.55556 

0.3535^ 
0.25140 
0.20298 


FIGURE  84 
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i  DELTA  T  =  1.2 

TFIN  =5.2 

T 

Q(N) 

1.2 
2.4 
3.6 

4.8 

0.38627 
0.49653 
0.27125 
0.21600 

FIGURE  85 
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FOR  DELTA  T  =  1.4 


T 


1.4 
2.8 
4.2 


TPIN  =5.2 

Q(N) 

0.41315 
0.42946 
0.22270 


FIGURE  86 
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FOR  DELTA 

T  = 

;  1. 

6 

TFIN  =5.2 

T 

Q(N) 

1.6 

3.2 

4.8 

0.41995 
O.36988 
0.19246 

FIGURE 

87 
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FOR  DELTA  T 

T 

1.8 
3.6 


=  1.8 


TFIN  =5.2 

Q(N) 

0.41369 

0.32118 


FIGURE  88 
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DELTA  T  =  2.0 

TFIN  =5.2 

T 

Q(N) 

2.0 
4.0 

0.40000 
0.28235 

FIGURE  89 
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D.   CASE  7  -  Y  +  |/tY  =  t2 

Following  the  procedure  as  In  II  C: 

f  ♦  FY  =  t2  P  =/t 


The  Laplace  transform  in  powers  of  s   : 


sy  +  Py  =  2/s3  y  =  — — — 


1  +  Ps"1 


Substituting  the  z-forms  of  Table  I  into  the  expression 
above  yields: 


yA*(z)  = 

T*1         z"1  +  kz   2  +  z    3 
6    '        (1-z-V 

720 

,    +   P   T        (l+z)"1 
(1  -  z) 

Dividing  by  T  and  rearranging  in  crescent  powers  of  z   : 

*(z)    =  -T3  +   (12^T3)z~1  +   (^T3)z~2  +  12*IT3z~3  =  T3z~lj  ,^-x 

A  360+l80FT-(l^0+360PT)z_1+2l60z~2-(1^0-360PT)z~3+ 

...  +  (360-l80PT)z_i< 


Choosing  for  first  value  for  T  =  0.2,  and  substituting 
in  the  previous  equation,  the  denominator  is  divided  into 
the  numerator,  after  the  value  of  P  is  replaced  in  the 
expression  above  for  the  corresponding  time. 
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The  results  and  the  plots  for  the  precise  and  approximated 
solution  follow,  with  iteration  times  from  0.2  sec.  to 
2.0  sec.  (Figs.  90  to  103).   Computer  algorithm  for  solution 
of  this  problem  is  in  Program  7  (page  128). 

DELTA  T  =  0.2  TFIN  =5.2 

T  Q(N)         PR 

0.0 
0.2 
0.4 
0.6 
0.8 
1.0 
1.2 
1.4 
1.6 
1.8 
2.0 
2.2 
2.4 
2.6 
2.8 
3.0 
3.2 

3.4 
3.6 
3.8 
4.0 
4.2 
4.4 
4.6 
4.8 
5.0 


■0.00002 

0.0 

0.00255 

0.00261 

0.01985 

0.02018 

0.06428 

0.06511 

0.14500 

0.14648 

0.26791 

0.27012 

0.43602 

0.43897 

0.65001 

0.65363 

0.90881 

0.91298 

1.21016 

1.21470 

1.55110 

1.55590 

1.92838 

1.93330 

2.33869 

2.34370 

2.77893 

2.78380 

3.24627 

3.25100 

3.73821 

3.74280 

4.25262 

4.25700 

4.78767 

4.79180 

5.34186 

5.34590 

5.91392 

5.91780 

6.50283 

6.50670 

7.10775 

7.11160 

7.72795 

7.73200 

8.36285 

8.36720 

9.01195 

9.01680 

9.67481 

9.68030 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C       Q(T)    =    T**2       &      P(T)    =SQRT(T) 
C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

INTEGER    *4ITB( 12)/12*0/ 
REAL    *4RTB(28)/28*0.0/ 

DIMENSION    A(IOO),     Q(IOO),    X(IOO),    XX(IOO),     PR(IOO) 
ITB(3)     =    5 
ITB(4)     =    5 
WRITE    (6,8) 
C 

DO    1    1=1,26 

READ    (5,5)     XX(I),PR(I) 

WRITE     (6,5)     XX(I),PR(I) 

1  CONTINUE 
C 

2  READ    (5,5,END=4)     TO.TFIN 
WRITE    (6,6)     TDtTFIN 
WRITE     (6,7) 

M    =    TFIN/TD 

T3    =    TD**3 

A(l)     =    -T3 

A(2)     =     124.0*T3 

A(3)    =    474.0*T3 

A(4)    =     124.0*T3 

A(5)     =    -T3 
C 

DO    3    N=1,M 

T    =    (N-1)*TD 

P    =    SQRT(T) 

PT    =    P*TD 

Fl    =    360. 0+180. 0*PT 

F2    =    -( 1440. 0+360. 0*PT) 

F3    =    2160.0 

F4    =    -( 1440. 0-360. 0*PT) 

F5    =    360. 0-180. 0*PT 
C 

Q(N)    =    A(N)/F1 
C 

A(N  +  1 )     =    A(N  +  1 )-(Q(N)*F2) 

A(N+2)     =    A(N+2)-( Q(N)*F3) 

A(N+3)     =    A(N+5)-(Q(N)*F4) 

A(N+4)     =    A(N+4)-(Q(N)*F5) 

AiN+5)     =    0.0 

WRITE     (6,5)     T,3(N) 

X(N)    =    T 

3  CONTINUE 
C 

ITB(i)     =    1 

ITB(2)     =    0 

ITBU2)     =    1 

CALL    DRAWP    (  26 , XX , PR , I TB  ,  RTB  ) 

ITB( 1)     =    3 

ITB(2)     =    2 

CALL    DRAWP    (  M  ,  X  ,  Q  ,  I  TB,  RTB  ) 

GO    TO    2 

4  STOP 
C 

5  FORMAT     (2F10.5) 

6  FORMAT     Cl't'FOR     DELTA    T=' , F6 . 2 , 4X , ■ TF I N= ' , F4. 1 ) 

7  FORMAT     ( •     T=« ,8Xf •     0(N)=' ) 

8  FORMAT     ( •        XX  =  '  , 7X,  •        PR='  ) 
END 


PROGRAM    7 
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x-scale  =  1.0  units/inch 
y-scale  =2.0  units/inch 

FOR  DELTA  T  =  0.2 


FIGURE  90 
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FOR  DELTA 

T  =  0 

3 

TFIN  =5.0 

T 

Q(N) 

0.0 

-0.00007 

0.3 

0.00832 

0.6 

O.O6328 

0.9 

0.19866 

1.2 

0.43235 

1.5 

0.76902 

1.8 

1.20437 

2.1 

1.72913 

2.4 

2.33229 

2.7 

3.00312 

3.0 

3.73225 

3.3 

4.51208 

3.6 

5.33673 

3.9 

6.20176 

4.2 

7.10386 

4.5 

8.04059 

FIGURE 

91 
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TE 


fin  i                                        r 

L  u  .1                                        »- 

22 

0G3                   oo; 

FOR  DELTA 

T 

=    0 

4 

TFIN      =5.0 

T 

Q(N) 

0.0 

-0.00018 

0.4 

0. 0189^ 

0.8 

0.14071 

1.2 

0.42726 

1.6 

0.89625 

2.0 

1.53631 

2.4 

2.32327 

2.8 

3.23134 

3.2 

4.23875 

3.6 

5.32927 

4.0 

6.49164 

4.1 

7.71831 

FIGURE 

92 
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FOR  DELTA 

T  = 

0.5 

TFIN  =5.0 

T 

Q(N) 

0.0 

-0.00035 

0.5 

0.035^1 

1.0 

0.25665 

1.5 

0.75346 

2.0 

1.52519 

2.5 

2.52824 

3.0 

3.71282 

3.5 

5.03967 

4.0 

6.48280 

4.5 

8.02638 

FIGURE 

93 
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FOR  DELTA 

T  = 

0. 

6    TFIN  =5.0 

T 

Q(N) 

0.0 

-0.00060 

0.6 

0.05842 

1.2 

0.41309 

1.8 

1.17315 

2.4 

2.29728 

3.0 

3.69936 

3.6 

5.30783 

4.2 

7.07854 

FIGURE 
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FOR  DELTA  T 

T 

0.0 
0.7 
l.ll 

2.1 

2.8 

3.5 
4.2 


=  0.7 


TFIN  =  5.0 

Q(N) 

-0.00095 
0.08844 
0.61018 
1.67849 
3.18952 
5.01324 
7.06634 


FIGURE  95 
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FOR  DELTA  T  =  0 . 8 


TFIN  =5.0 
Q(N) 


0.0 
0.8 
1.6 
2.4 
3.2 
4.0 


-0. 
0, 
0, 
2 
4 
6 


00142 
12570 
84689 
26036 
18196 
44458 


FIGURE  96 
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DELTA  T   =   0.9 

TFIN  =5.0 

T 

Q(N) 

0.0 

0.9 
1.8 

2.7 
3.6 

-0.00202 
0.17030 
1.12150 
2.90981 
5.25900 

FIGURE    97 
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FOR  DELTA  T  =  1.0 

T 

0.0 
1.0 
2.0 
3.0 
4.0 

FIGURE  9  8 


TFIN  =5.0 

Q(N) 

-0.00278 
0.22222 
1.43192 
3.61877 
6.40907 
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FOR  DELTA  T 

T 

0.0 
1.2 

2.4 
3.6 


=  1.2 


TFIN  =5.0 

Q(N) 

-0.00480 
0.34756 
2.15138 
5.18870 


FIGURE  99 
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FOR  DELTA 

T  = 

1. 

i] 

TFIN  =5.0 

T 

Q(N) 

0.0 
l.ll 

2.8 

-0.00762 
0.50030 
2.98831 

FIGURE 

100 
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FOR  DELTA  T  =  1.6 

T 

0.0 
1.6 

3.2 


TFIN  =5.0 

Q(N) 

-0.01138 
0.67862 
3.92797 


FIGURE  101 
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FOR  DELTA  T   =    1.8        TFTN   =    5.0 
T 


0.0 
1.8 


Q(N) 

-0.01620 
0.88064 


FIGURE   102 
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FOR  DELTA  T 

T 

0.0 
2.0 


=  2.0 


TFIN  =5.0 

Q(N) 

-0.02222 
1.10457 


FIGURE  10  3 
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1.      Case    8  -  Y  +   tY  =   t2 

Same    calculations   as   in   case   7  with  exception   for 
the   substitution   on   the  value   of  P   during  the   division 
steps,    of  Equation    (47),    repeated  here,   which  now  has   the 
value   P  =   t . 

The   algorithm  for   computer   solution    (Program   8, 
page    144) ,    the   output    for   several   time   iterations   and  the 
respective   graphs    follow    (Figs.    104   to   117) 


*,v  =       -T3+124T3)z~1+(474T3z~2+124T3z~3-T3z~It 

A  360+l80PT-(l440+360PT)z_1+2l60z~2-(l440-360PT)z"3+. . . 

-4 


FOR   DELTA   T   =    0.2 


0.0 
0.2 
0.4 
0.6 
0.8 
1.0 
1.2 
1.4 
1.6 
1.8 
2.0 
2.2 
2.4 
2.6 
2.8 
3.0 
3.2 
3.4 
3.6 
3.8 
4.0 
4.2 
4.4 
4.6 
4.8 


•   • 

•+(360-l80PT)z' 

TFIN   =5.0 

Q(N-) 

PR 

0.00002 

0.0 

0.00261 

0.00264 

0.02041 

0.02066 

0.06629 

0.06707 

0.14905 

0.15068 

0.27254 

0.27522 

0.43567 

0.43942 

0.63339 

O.638OI 

0.85809 

0.86328 

1.10124 

1.10660 

1.35473 

1.35990 

1.61184 

1.61660 

I.86765 

1.87190 

2.11908 

2.12280 

2. 36463 

2.36780 

2.60393 

2.60670 

2.83735 

2.83980 

3.06563 

3.06780 

3.28958 

3.29160 

3.51003 

3.5H90 

3.72765 

3.72940 

3.94301 

3.94460 

4.15655 

4.15810 

4.36861 

4. 37000 

4.57945 

4. 58080 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C      Q(T)    =    T**2       £      P(T)     =T 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

INTEGER    *4ITB( 12)/12*0/ 
REAL    *4RTB(28)/28*0.0/ 

DIMENSION  A(IOO),  Q(IOO),  X(IOO),  XX(IOO),  PR(IOO) 
ITB(3)  =  5 
ITB(4)  =  5 
WRITE  (6,8) 
C 

DO    1    1=1,26 

READ    (5,5)    XX( I) ,PR(I ) 

WRITE     (6,5)     XX(I),PR(I) 

1  CONTINUE 
C 

2  READ    (5,5,END=4)     TD.TFIN 
WRITE     (6,6)     TD,TFIN 
WRITE    (6,7) 

M    =    TFIN/TD 

T3    =    TD**3 

A(l)    =    -T3 

A(2)     =     124.0*T3 

A(3)    =    474.0*T3 

A(4)     =     124.0*T3 

A(5)    =    -T3 
C 

DO    3    N=1,M 

T    =    (N-1)*TD 

P    =    T 

PT    =    P*TD 

Fl    =    360. 0+180. 0*PT 

F2    =    -( 1440. 0+360. 0*PT) 

F3    =    2160.0 

FA    =    -( 1440. 0-360. 0*PT) 

F5    =    360. 0-180. 0*PT 
C 

Q(N)    =    A(N)/F1 
C 

A(N+1 )     =    A(N+1)-(Q(N)*F2) 

A(N+2)     =    A(N+2)-( Q(N)*F3) 

A(N+3)     =    A(N+3)-(Q(N)*F4) 

A(N+4)     =    A(N+4)-(Q<N)*F5) 

A(N+5)     =    0.0 

WRITE    (6,5)     T,Q(N) 

X(N)    =    T 

3  CONTINUE 
C 

ITB(l)     =    1 

ITB(2)     =    0 

ITB(12)     =    1 

CALL    DRAWP    ( 26 , XX , PR , I Tb , RTB ) 

ITB( 1)     =    3 

ITB(2)     =    2 

CALL    DRAWP    (  M , X , Q , I TB , RTB ) 

GO  TO  2 

4  STOP 
C 

5  FORMAT  (2F10.5) 

6  FORMAT  (•lS'FOR  DELTA  T=  '  ,  F6  .  2  ,  4X  ,  '  TF  I  N  =  '  ,  F4  . 1  ) 

7  FORMAT  ( •  T=' ,8X, '  0(N)=' ) 

8  FORMAT 
END 


PROGRAM  8 
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x-scale  =  1.0  units/inch 
y-scale  =1.0  units/inch 


FOR  DELTA  T 


0.2 
FIGURE  104 
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FOR  DELTA 

T  =  0. 

3 

TFIN  =5.0 

T 

Q(N) 

0.0 

-0.00007 

0.3 

0.00861 

0.6 

0.06534 

0.9 

0.20305 

1.2 

0.43105 

1.5 

0.73670 

1.8 

1.09444 

2.1 

1.47676 

2.4 

1.86218 

2.7 

2.23828 

3.0 

2.60051 

3-3 

2.94932 

3.6 

3.28723 

3.9 

3.61706 

4.2 

3.94109 

4.5 

4.26094 

FIGURE 

105 
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DELTA  T  =  0 

.4 

TFIN  =5.0 

T 

Q(N) 

0.0 

-0.00018 

0.4 

0.01975 

0.8 

0.14440 

1.2 

0.42470 

1.6 

0.84250 

2.0 

1.33874 

2.4 

1.85444 

2.8 

2.35490 

3.2 

2.83016 

3.6 

3.28380 

4.0 

3.72261 

4.4 

4.15198 

FIGURE 

106 
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FOR  DELTA 

T  =  0 

.5 

TFIN  =5.0 

T 

Q(N) 

0.0 

-0.00035 

0.5 

0.03704 

1.0 

0.25926 

1.5 

0.71717 

2.0 

1.32660 

2.5 

1.97228 

3.0 

2.58929 

3.5 

3.16744 

4.0 

3.71877 

4.5 

4.25487 

FIGURE 

107 
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FOR  DELTA  T  =  0.6     TFIN  =5.0 


T 


0.0 
0.6 
1.2 

1.8 
2.4 
3.0 
3.6 
4.2 


Q(N) 

-0.00060 
0.06102 
0.40738 
1.05761 
1.83168 

2.58151 
3.27410 
3.93011 


FIGURE  108 
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FOR  DELTA 

T  =  0.7 

TFIN  =5.0 

T 

Q(N) 

0.0 

-0.00095 

0.7 

0.09183 

1.4 

O.58367 

2.1 

1.42363 

2.8 

2.32707 

3.5 

3.15544 

4.2 

3.92482 

FIGURE 

109 
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FOR   DELTA   T   -    0.8        TFIN    =5.0 
T  Q(N) 


0.0 

-0.00142 

0.8 

0.12929 

1.6 

0.78206 

2.4 

1.79807 

3.2 

2.80113 

4.0 

3.70243 

FIGURE  110 
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FOR   DELTA   T   =    0.9        TFIN   =5.0 
T  Q(N) 


0.0 

=0.00202 

0.9 

1.8 

2.7 
3.6 

0.17295 
0.99663 
2.16991 
3.25361 

FIGURE  111 
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FOR  DELTA  T  =  1.0 

T 

0.0 
1.0 
2.0 
3.0 
4.0 

FIGURE  112 


TFIN  =5.0 

Q(N) 

-0.00278 
0.22222 
1.22222 

2.53333 
3.68888 
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Dl 

SLTJ 

\   T   =    1., 

2 

TPIN  =    5. 

0 

T 

Q(N) 

0 
1 

2 

3. 

.0 
.2 

6 

FIGURE 

113 

-0.00480 
0.33488 
1.69089 
3.22784 
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FOR  DELTA  T  - 

1 

4 

TFIN  =5.0 

T 

Q(N) 

0.0 
1.4 

2.8 

-0.00762 

0.46195 
2.16618 

FIGURE 

114 
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1  DELTA  T  =  1.6 

TPIN  =5.0 

T 

Q(N) 

0.0 
1.6 

3.2 

-0.01138 
0.59883 
2.63754 

FIGURE  115 
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FOR  DELTA  T  =  1.8    TFIN  =5.0 
T 


0.0 
1.8 


Q(N) 

-0.01620 
0.74198 


FIGURE  116 
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1 

::z 

aa 

oai 

FOR  DELTA   T   =    2. 

0 

TFIN   = 

5. 

0 

T 

Q(N) 

0.0 
2.0 

FIGURE 

117 

-0.02222 
0.88889 
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2.   Case  9  -  Y  +  t2Y  =  t2 

Same  reasoning  as  in  case  8.   This  time  P  =  t2. 

Computer  algorithm  for  solution  of  this  problem  is 
shown  in  Program  9  (Page  160) .   Solutions  and  graphs  for 
several  T's  follow  in  Figures  118  to  131. 

FOR   DELTA   T   =    0.2 


0.0 
0.2 
0.4 
0.6 
0.8 
1.0 
1.2 
1.4 
1.6 
1.8 
2.0 
2.2 
2.4 
2.6 
2.8 
3.0 
3.2 
3.4 
3.6 
3.8 
4.0 
4.2 
4.4 
4.6 
4.8 


TFIN 

=  5.0 

Q(N) 

PR 

■0.00002 

0.0 

0.00266 

0.00266 

0.02098 

0.02110 

0.06883 

0.06947 

0.15509 

0.15690 

0.27985 

0.28347 

0.43228 

0.43785 

0.59256 

0.59934 

0.73812 

0.74469 

0.85183 

O.85686 

0.92749 

0.93051 

0.96977 

0.97125 

0.98924 

0.99002 

0.99645 

0.99714 

0.99857 

0.99934 

0.99909 

0.99988 

0.99924 

0.99998 

0.99932 

1.00000 

0.99938 

1.00000 

0.99943 

1.00000 

0.99947 

1.00000 

0.99950 

1.00000 

0.99952 

1.00000 

0.99954 

1.00000 

0.99956 

1.00000 

(z) 


Equation    (47)    is    repeated   for   convenience 
-T3+(124T3)z-1+(WT3)z~2+124T3z"3-T3z"1, 


A  360+l80PT-(l440+360PT)z~1+2l60z  2-(l44o-360PT)z  3+- 

...  +  ( 360-1 80PT)z-ij 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C       Q(T)    =    T**2      £      P(T)    =  T**2 
C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

INTEGER    *4ITB( 12)/12*0/ 
REAL    *4RTB(28)/28*0.0/ 

DIMENSION  A(IOO),  Q(IOO),  XUOO),  XX(IOO),  PR(IOO) 
ITB<3)  =  5 
ITB(4)  =  5 
WRITE  (6,8) 
C 

DO    1    1=1,26 

READ    (5,5)    XX( I ) ,PR( I ) 

WRITE    (6,5)    XX( I) ,PR( I) 

1  CONTINUE 
C 

2  READ    (5,5,END=4)     TD,TFIN 
WRITE    (6,6)     TD,TFIN 
WRITE    (6,7) 

M    =    TFIN/TD 

T3    =    TD**3 

A(l)    =    -T3 

A(2)    =     124.0*T3 

A(3)    =    474.0*T3 

A(4)    =     124.0*T3 

A(5)    =    -T3 
C 

DO    3    N=1,M 

T    =    (N-1)*TD 

P    =    T**2 

PT    =    P*TD 

Fl    =    360. 0+180. 0*PT 

F2    =    -< 1440. 0+360. 0*PT) 

F3    =    2160.0 

F4    =    -( 1440.C-360.0*PT) 

F5    =    360. 0-180. 0*PT 
C 

Q(N)    =    A(N)/F1 
C 

A(N+1)     =    A(N+1)-(Q(N)*F2) 

A(N+2)     =    A(N+2)-(Q(N)*F3) 

A(N+3)     =    A(N+3)-( Q(N)*F4) 

A(N+4)     =    A(N+4)-(Q(N)*F5) 

A(N+5)     =    0.0 

WRITE     (6,5)     T,Q(NJ 

X(N)    =    T 

3  CONTINUE 
C 

ITB(l)     =    1 

ITB(2)     =    0 

ITB(12)     =    1 

CALL    DRAWP    ( 26 , XX , PR ,  I  TB, RT B ) 

ITB( 1)     =3 

ITB(2)     =    2 

CALL    DRAWP    ( M , X , Q , I TB , RTB ) 

GO  TO  2 

4  STOP 
C 

5  FORMAT  (2F10.5) 

6  FORMAT  CI', 'FOR  DELTA  T=',F6.2 

7  FORMAT  ( '  T=« ,8X,  '  Q(N)  =  '  ) 

8  FORMAT  ( '   XX=' ,7X,'   PR=«) 
END 


PROGRAM  9 
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III! — h-H — ^  +  +  -I — h-t- 


x-scale  =1.0  units/inch 
y-scale  =0.2  units/inch 

FOR  DELTA  T  =  0.2 


FIGURE  118 
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FOR  DELTA   T   =    0.3 

TFIN   =5.0 

T 

Q(N) 
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0.6 
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1.2 

0.42549 

1.5 

0.65988 

1.8 

0.84534 

2.1 

0.94944 

2.4 

0.98840 

2.7 

0.99809 

3.0 

0.99819 

3.3 

0.99851 

3.6 

0.99877 

3.9 

0.99895 

4.2 

0.99910 

4.5 

0.99921 

FIGURE 
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FOR   DELTA 

T    = 

0.4 

TFIN    =5.0 

T 

Q(N) 

0.0 

-0.00018 

0.4 

0.02067 

0.8 

0.15013 

1.2 

0.41634 

1.6 

0.71810 

2.0 

0.91765 

2.4 

0.98739 

2.8 

0.99659 

3.2 

0.99714 

3.6 

0.99787 

4.0 

0.99827 

4.4 

0.99859 

FIGURE    120 
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FOR  DELTA 

T    = 

0. 

5 

TFIN   =5.0 

T 

Q(N) 

0.0 

-0.00035 

0.5 

0.03922 

1.0 

0.26275 

1.5 

0.63278 

2.0 

0.90925 

2.5 

0.99187 

3.0 

0.99499 

3-5 

0.99641 

4.0 

0.99731 

4.5 

0.99789 

FIGURE 
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FOR   DELTA  T   = 

0. 

6        TFIN   =5.0 

T 

Q(N) 

0.0 

-0.00060 

0.6 

0.06498 

1.2 

0.39243 

1.8 

0.80675 

2.4 

0.98482 

3.0 

0.99326 

3.6 

0.99498 

4.2 

0.99658 

FIGURE 

122 
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FOR  DELTA 

T  = 

0. 

7   TFIN  =5.0 

T 

Q(N) 

0.0 

-0.00095 

0.7 

0.09760 

1.4 

0.52265 

2.1 

0.91859 

2.8 

0.99655 

3.5 

0.99033 

4.2 

0.99646 

FIGURE 
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FOR   DELTA 

T   =    C 

).8 

TFIN  =5.0 

T 

Q(N) 

0.0 
0.8 
1.6 
2.4 
3.2 
4.0 

-0.00142 
0.13588 
0.64020 

0.97679 
0.98920 
0.99299 

FIGURE 
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FOR  DELTA  T  =  0.9 


TFIN  =5.0 


0.0 
0.9 
1.8 

2.7 
3.6 


Q(N) 

-0.00202 
0.17809 
0.73807 
0. 99964 
0.9823^ 


FIGURE  125 
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FOR   DELTA 

T   =    1 

.0 

TFIN   =5.0 

T 

Q(N) 

0.0 
1.0 
2.0 
3.0 

M.O 

-0.00278 
0.22222 
0.81481 
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FIGURE 
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FOR  DELTA   T 

T 

0.0 
1.2 

2.4 
3.6 


=    1.2 


TFIN    =5.0 

Q(N) 

-0.00480 
0.30901 
0.91^28 
0.99H7 


FIGURE   127 
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FOR  DELTA  T  =  1.4 


0.0 

1.4 

2.8 


TFIN  =5.0 

Q(N) 

-0.00762 

0. 38561 
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FIGURE   128 
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FOR   DELTA 

T   = 

1. 

6           TFIN   =5.0 

T 

Q(N) 

0.0 
1.6 
3.2 

-0.01138 
0.  -144794 

0.98867 

FIGURE 
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FOR  DELTA   T 

T 

0.0 
1.8 


=   1.8 


TFIN  =5.0 

Q(N) 

■0.01620 
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FIGURE    130 
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FOR  DELTA  T  =  2.0    TFIN  =  5.0 
T  Q(N) 


0.0 
2.0 


-0, 
0, 


02222 
53333 


FIGURE  131 
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E.   CASE  10 

No  one  of  the  previous  cases  has  a  time  varying 
coefficient  which  is  exponential.   So,  from  Ref.  2  this 
case  was  extracted.   The  equation  is: 


C  ||  +  ^^7  =  0       being  r(t)  =  R  exp  (^-) 


where  T,  is  a  factor  in  the  variation  of  r.   Dividing  all 
the  terms  by  C  yields: 


||  +  -L  exp(^.>  e  =  0 


where  T2  =  RC . 

For  this  example  T,  and  T~  were  made  equal  to  one  and 
6/._n>>  =  1.0  .  Rewriting  the  equation  and  expressing  its 
Laplace  transform  in  powers  of  s~  : 


e  +  Ce  =  0        C  =  et 


sE  -  1  +  CE  =  0 

E    =        1       ~    1  CS"2 


S  +  C    s    1  +  Cs-1 


If  the  rules  in  Example  2,  Section  IIB  were  followed, 
substituting  the  z-forms  of  s~k,  dividing  by  T  and  adding 
the  final  results  of  the  quotient  to  1.0  which  is  the  time 
function  for  1/s  : 
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r   li  1 + lOz"1 +  z~2 
...    12   (l-z-1)2 
A    i  +  c  T  1  +  z"1 

1   0?l-2-l 


For  T  =  0.1 

-1        -2 
E  *  _  -0. 1  -  z  x  -  O.lz 


A     (12/C  +  0.6)  -  (2VOz_1  +  (12/C  -  0.6)z"2 

Dividing  the  denominator  into  the  numerator  and  adding  one 
to  the  result  one  gets 


0.99206  +  0.91i*12z  1   +  0.85007z  2  + 


However,  looking  at  the  graph  (Fig.  132)  the  precise 
and  the  approximated  solution  start  diverging  after  a  while 
This  is  due  to  the  fact  that  in 


s  +  C 
if  C  goes  to  zero  the  equation  becomes 

E  =  1/s 

which  is  the  final  value  that  the  approximated  solution 
approaches. 
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This  fact  is  very  important  because  it  leads  to  these 
important  conclusions: 

1.  If  C.  goes  to  zero  with  time,  first  check  the  Laplace 
transform  and  the  initial  equation  substituting  C^  by  zero. 
If  by  inspection  they  do  not  lead  to  the  same  type  of 
solution,  this  method  does  not  apply. 

2.  Verify,  before  doing  the  division  if  any  of  the 
varying  coefficients,  unexpectedly,  goes  to  infinity;  by 
dividing  by  zero,  for  instance. 


H 

o 


+  +  +  + 


-+-  +  +  +  +  +  +  + 


oo ; 


OJ'o 


Figure  132 


X-scale 
Y-scale 


1.0  units/inch 
.5  units/inch 


—  Precise  solution 
+  Approximate  solution 
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IV.   CONCLUSIONS 

This  technique,  for  the  successfully  area  covered  in 
this  paper  (linear  differential  equations  with  time  varying 
coefficients),  showed  to  have  a  good  precision  even  for 
relatively  large  iteration  steps.   The  labor  necessary  to 
complete  a  problem  by  handwork  is  not  very  involved  and 
needs  only  a  small  algorithm  for  computer  solution  with 
C.P.U.  times  much  smaller  than  existing  numerical  methods. 

However,  with  time  varying  coefficients  that  go  to  zero 
with  time,  attention  must  be  paid  to  the  application  of 
this  method.   As  a  general  rule,  this  method  does  not  apply. 

In  some  cases,  the  final  equation  for  the  approximate 
solution,  before  the  division  is  performed,  must  be  altered 
in  order  to  not  have  time  varying  coefficients  (referred  as 
C.  in  Section  IIC)  in  the  numerator.   Caution  is  necessary 
to  verify,  before  starting  the  division  steps  if  any  varying 
coefficient  unexpectedly  has  a  value  of  infinity;  by  division 
by  zero  when  substituting  a  particular  value  of  C^. 

The  technique  tested  in  this  paper  does  not  apply  to  the 
case  of  linear  differential  equations  with  varying  coefficients 
only.   Linear  differential  equations  with  constant  coefficients 
and  nonlinear  differential  equations  are,  as  well,  solved  by 
this  technique. 

These  areas,  in  the  opinion  of  the  author,  are  valuable 
as  topics  for  other  theses. 
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